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ABSTRACT 

We present the results of precision gamma-ray timing measurements of the binary millisecond pulsar 
PSR J2339-0533, an irradiating system of “redback” type, using data from the Fermi Large Area Telescope. 
We describe an optimized analysis method to determine a long-term phase-coherent timing solution spanning 
more than six years, including a measured eccentricity of the binary orbit and constraints on the proper motion 
of the system. A major result of this timing analysis is the discovery of an extreme variation of the nom¬ 
inal 4.6 hr orbital period Path over time, showing alternating epochs of decrease and increase. We inferred 
a cyclic modulation of Porb with an approximate cycle duration of 4.2 years and a modulation amplitude of 
APorb/ Porb = 2.3 X 10 Considering different possible physical causes, the observed orbital-period modula¬ 
tion most likely results from a variable gravitational quadrupole moment of the companion star due to cyclic 
magnetic activity in its convective zone. 

Subject headings: gamma rays: stars - methods: data analysis - pulsars: individual (PSR J2339-0533) 


1. INTRODUCTION 


The Large Area Telescope (EAT; [Atwood et al.||2009) on 
the Fermi satellite has paved the way to greatly increase the 
known Galactic population of compact binary systems har¬ 
boring irradiating millisecond pulsars (MSP). These intrigu¬ 
ing systems are commonly re ferred to as “black widows” 
and “redbacks” ( [Roberts [2013[ ), as the intense pulsar radia¬ 
tion is gradually destroying the companion star. Typically, 
in black widows the mass of the degenerate companion is 
very low (~ 0.008-0.05M©), whereas the redback type is 
distinguished by having a non-degenerate, main-sequence- 
like companion that is more massive (^ O.lS-O.bM©). Both 
types of systems provide interesting opportunities to study the 
interaction between the pulsar wind and the stellar compan¬ 
ion, their unusu al formation history ( [Chen et al.|2013[ Smed- 
[ley et al.][2015[ ) an d the evolutionary link to low-mass X-ray 
binaries (LMXBs) ( Archibald et al. [20091 Papitto et al. 2013|l, 
as well as the masses of neutron kars ^van Kerkwijk et al.[ 
[201 l[[Romani et al.|2012| l. 

Since the launch of Fermi, numerous new black widow and 
redback pulsars were found in targeted pulsar searches of for¬ 
merly unidentified EAT gamma-ray sources with radio tele- 


scopes (e.g., [Ransom et al.|2011[ Ray et al.|2012[ Barr et al. 
[2013|l. Only in one case a direct search in LAT data revealed 


the ga mma-ray pulsations of PSR J1311-3430 ( [Pletsch et al. 
[2012J by exploiting partial knowledge of the orbit obtained 
from prior identification of the heated c ompanion at optical 
and X-ray wavelengths (|Romani[[2012[l. Although the pul¬ 
sar was subsequently detected in the radio band ([Ray et al. 
[2013| l, it was shown to be weak and intermittent, which had 
precluded a detection in typical radio searches. 

The discovery path to the redback system PSR 12339-0533 
followed a similar route. The formerly unidentified LAT 
source OFGL J2339.8-0530 was one of the bright gamma- 
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ray sour ces unveiled durin g the first 3 months of the Fermi 
mission ( Abdo et al.|2009| l. X-ray and optical follow-up ob- 
servations identihed a probable counterpart, which showed 
many properties reminiscent of an MSP in a 4.6 hr binary 
orbit around a low-mass companion ([Romani & Shaw|2011| 
Kong et al.|2012[|. Only pulsations were lacking for an unam¬ 
biguous conhrm ation of the pulsar natur e. Using the Green 
Bank Telescope, [Ray et al.[([2014[[2015|l detected 2.8 ms ra¬ 
dio pulsations. They showed thait the companion was sub¬ 
stantially more massive than modeling of the optical data 
had initially suggested. With the spin period and orbital pa¬ 
rameters tightly constrained, they were also able to detect 
gamma-ray pulsations with the LAT, verifying the identifica¬ 
tion of OFGL J2339.8-0530. Their observed dispersion mea¬ 
sure (DM) from the rad io detection provided an estimated dis¬ 
tance of about 450pc ( Ray et ar][2014 2015| l. Up to now, a 
conjectured variability in the orbital parameters had precluded 
a phase-coherent ephemeris for the pulsar covering the entire 
time span of the available LAT data. 

Other black widow and redback MSP systems have also 
been obser ved to show significant variations in their orbital 
parameters ([Arzoumanian et al.|1994][Doroshenko et al.[2001[ 


Archibald et al.|2013[ l. In the interest of better understanding 

such phenomena, it is essential t o carry out precis i on timing 
over longer time intervals (e.g., [Nice et al.[|2000[ [Lazaridis| 
et al. 201 1[ ). Often the eclipse of the pulsar’s radio emis- 
sion over a large fraction of the orbit additionally complicates 
or even prohibits an accurate timing analysis of the orbital- 
parameter variability. On the contrary, the pulsar’s gamma- 
ray emission is essentially unaffected by this problem. Thus, 
the continuously recorded multi-year LAT survey-mode data 
is uniquely suited to carefully monitor the long-term evolution 
of the orbital parameters of those systems. 

Here, we present the results of a precision timing analysis 
of PSR J2339-0533 using LAT gamma-ray photon data that 
cover more than six years. In Section!^ we describe the data 
preparation, the improved timing methodology, and provide 
the details of the inferred phase-coherent pulsar ephemeris. 
Section [^provides a detailed discussion of these results. Fi¬ 
nally, a summarizing conclusion follows in Section]^ 
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2. GAMMA-RAY DATA ANALYSIS 
2.1. Data Preparation 

For this analysis we employed Pass 7 reprocessed LAT pho¬ 
ton data recorded between 2008 August 4 and 2014 Octo¬ 
ber 17. Via the Fenni Science Tool^we selected photons be¬ 
longing to the P7REP_SOURCE class with reconstructed di¬ 
rections within 8° of the pulsar position, energies from 0.1 to 
100 GeV, and zenith angles < 100°. Photons recorded when 
the LAT’s rocking angle exceeded 52° or when the LAT was 
not in nominal science mode were excluded. 

We constructed spectral models including all gamma-ray 
sources found within 13° of the nominal pulsar position from 
a preli minary version of th e third Fermi LAT source catalog 
(3FGL; Acero et al.|2015| l based on four years of LAT data. 
The parameters were left free only for point sources within 
the inner 5°. The source models included contributions from 
the Galactic diffuse emission, the extragalactic diffuse emis¬ 
sion, and the residual instrumental backgrouncjj The pulsar 
spectrum was modeled as an exponentially cutoff power law, 
dN/dE oc E~^ exp (-£ /Ec ), where F denotes the spectral in¬ 
dex and Ec is the cutoff energy. Using this spectral model 
of the region and the P7REP_SOURCE_V15 instrument re¬ 
sponse functions, we em ployed gts rcprob to compute a 
weight for each photon (|Kerr||201 l|l, measuring the proba¬ 
bility of having originated from the pulsar for ftirther back¬ 
ground suppression. 

From the final data set we discarded photons with probabil¬ 
ity weights less than 0.001. While this chosen weight cutoff 
makes only a negligible change to the significance of the de¬ 
tected pulsations, it dramatically improves the computational 
efficiency for the subsequent timing analysis. 

2.2. Timing Analysis 

A central aspect of this work is to obtain a precise 
ephemeris, i.e. a phase-coherent timing solution, covering the 
full time range of LAT data available. This means counting 
the exact number of pulsar turns over the past six years of 
the Fermi mission. The initial ephemeris obtained from the 
radio pulsar discovery using a circular-orbit mode l was suf- 
ficient to reveal significant gamma-ray pulsations (|Ray et al. 
2014|l, but was unable to maintain accurate phase coherence 


over a time span longer than about three years, most likely 
due to unaccounted orbital variability of the PSR J2339-0533 
system. This motivated the present investigation to study and 
model the putative orbital parameter variation, facilitating a 
long-term phase-coherent timing solution. 

To begin with, we examined the variation of the orbital pe¬ 
riod Porb and the projected semimajor axis x across neigh¬ 
boring subsets of data. For this, we fixed the sky position 
to the location of the optical counterpart at a = 23*'39™38?75 
and 5 = -05°33'05"3 (J2000.0) from |Romani & Shaw| ( [Ml 


The remaining pulsar parameters were set to those of the pre¬ 
liminary ephemeris. Assuming a circular orbit we scanned 
ranges in {PoibA} on a dense grid around their initial values 
at fixed time of ascending node Tasc- At each grid point we 
computed the weighted //-test statistic ( |de Jager et al.|[T98^ 
|Kerr|20lT] l using photons within a fixed time window of size 
110 days. This time window covered almost two preces¬ 
sion periods of Fermi (56 days) and was chosen by balanc¬ 
ing signal-to-noise ratio and time resolution, being just long 

http://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/overview.html 
^ http://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html 


enough to still accumulate a detectable signal-to-noise ratio. 
This window was then slid over the entire data set using a 
50% overlap between subsequent steps. 

The results of this study provided the first evidence that 
PSR J2339-0533 undergoes alternating intervals of orbital- 
period increase and decrease, whereas for x no significant 
variation was observed. To obtain an approximate solution for 
the evolution of Porb over time, we considered the Porb value 
giving the highest //-test in each sliding-window step. It be¬ 
came apparent that modeling this set of Porb values over time 
required us to employ a polynomial expansion of the orbital 
frequency /orb = 1 /Porb about Pasc^ so that for a given time t 
the orbital phase /orb (measured in cycles) is written as 


Dib(0 — ^ ^ 


k=a 


Ak) 

J orb 

(k+l)! 


(f-TAsc)'"^'^ , 


( 1 ) 


where represents the kth time-derivative of /orb. Simple 
least-square minimization suggested that the order K needed 
to give an adequate fit was at least five or six. Ultimately, to 
determine the model providing the best fit we used the more 
sensitive timing analysis described next. 

As a starting point for the precision timing analysis, we ex¬ 
tended the initial pulsar ephemeris by including six orbital- 
frequency derivatives based on the approximation for the 
orbital-period evolution of PSR J2339-0533 obtained in the 
previous step. Since the orbit is apparently of extremely low 
ecce ntricity e, we used the Lagrange-Laplace parametriza- 
tion ( |Lange et al. 200 1| | to describe the binary motion with 
ei = esinw and t 2 = ecosw, such that e = + The 

co nstruction of a pul se profile template foll owed the method 
of |Abdo et al. |( |2013] l extending the work of |Ray et al.| ( l201 l| l 
by also including the photon weights, labeled Wj for the jth 
photon. Since this requires a valid ephemeris, we first used 
only the subinterval (about half) of the full data set over 
which the initial ephemeris maintained phase coherence for 
this purpose. In general the pulse profile F{cj)) can be rep¬ 
resented by a wrapped probability density function (PDF) of 
the pulsar’s rotational phase / G [0,1), measured in turns. For 
PSR J2339-0533, we found that this PDF is well approxi¬ 
mated by a sum of a constant bac kground and five u nimodal 
wrapped Gaussian distributions (cf.|Abdo et al.|2013||. The ro¬ 
tational phase /(f/, u) is a function of the photon arrival time tj 
and the vector u that collects all pulsar spin, positional, and 
orbital parameters. Thus, to obtain the best-fit values f or u of 
a given timing model, we evaluated the log-likelihood (lAbdo' 
let al.|2013| l, 

N 

\og£{u) = '^log[wjFi<j)j{tj,u)) + (l-Wj)] , (2) 

/=! 

where N is the total number of gamma-ray photons. 

To find the parameters that give the global maximum of the 
log-likelihoo d, we followed a n ovel approac h that differs from 
the method ofjRay et al.|(j201 l|l. Note that in|Ray et al.|P011 1 
the LAT data were subdivided into segments in time and over 
each such segment the photon times were folded using a pre¬ 
liminary ephemeris to obtain a pulse time-of-arrival (TOA) 
measurement. The set of TOAs from all segments then served 
as input to a global fitting procedure. Due to the sparseness 
of the gamma-ray photons, a sufficiently high signal-to-noise 
ratio for a single TOA determination required long integration 
times of order weeks. However, this makes it very difficult to 
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Exposure Rotational phase Residual (//s) AP^^bZ-Porb {^0”^) 

(10® cm® s) 


Fig. 1.— Panel A: exposure versus time, binned in steps of 14 days. Panel B: integrated weighted pulse profile, using 50 bins per rotation. The vertical 
error bars show statistical Icr uncertainties. Two rotational cycles are shown for clarity. Panel C: two-dimensional weighted histogram of rotational phase versus 
time, using 140 bins over time. The weighted counts are not exposure-connected. Panel D: formal timing residuals, where the error bars show statistical Icr 
uncertainties. The data points were obtained from non-overlapping subsegments of data that gave approximately the same signal-to-noise ratio (i.e. an P-test 
value of about 15). Only significant measurements with a log-likelihood greater than 8 were included. Panel E: time evolution of the orbital-period change (solid 
black curve) based on the timing solution of Tableflland the surrounding gray-shaded region shows the statistical Icr uncertainties. The uncertainties increase with 
the time distance to Tasc (shown by the horizontaTaotted line), since the polynomial expan sion of the orbital-frequency used in the timing model is about Tasc- 
The dashed curve represents the best-fit cyclic modulation model discussed in Section [T^ and the dotted-dashed line shows the linearly changing component of 
this model. 


detect effects on much shorter time scales, and orbital param¬ 
eter variability in particular, since the preliminary ephemeris 
used to fold the photon times in each segment is not necessar¬ 
ily accurate. Given the pulsed flux level of PSR J2339-0533, 
one would need segments of duration at least 30 days. For 
comparison, this interval spans already about 155 orbital rev¬ 
olutions. 

Therefore, we instead evaluated the log-likelihood of Equa¬ 
tion 0 directly using the entire set of unfolded gamma-ray 
photon times (without binning in time). For the exploration of 
log£ over the relevant parameter space of the timing model 


we employed the Monte Carlo (MC) sampler MultiNest 
( Feroz et al.||2009 1 . This multimodal nested sampling algo- 
rithm is especially efficient in sampling challenging likeli¬ 
hood surfaces and allows one to calculate posterior distribu¬ 
tions as a byproduct. 

The timing analysis remains an iterative procedure, specifi¬ 
cally since the actual pulse profile is not known a priori. In the 
first MC run we still held the sky position fix ed to the optical- 
counterpart location (|Romani & Shaw|20111l, but let the spin 
and orbital parameters vary. ITiis gave a hrst phase-connected 
solution spanning the full FAT data set, from which we then 
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TABLE 1 

Parameters for PSR J2339-0533 


Parameter Value 

Range of observational data (MJD). 54683 - 56947 

Reference epoch (MJD). 55100.0 


Timing Parameters 

R.A., a (J2000.0). 23^39"'38?74(1) 

Decl., 5 (J2000.0). -05°33^05f'32(3) 

Proper motion in a,/io; cos <5 (mas yr“^). 11(4) 

Proper motion in 5, fis (mas yr“^). -29(10) 

Spin frequency, / (Hz). 346.71337922051(2) 

1st spin frequency derivative, / (Hz s“^). —1.6952(8) x 10“^^ 

2nd spin frequency derivative, |/| (Hz s“^) .... < 2 x 10“^^ 

Orbital frequency, /o^b (10^^ Hz). 5.993873572(1) 

1st derivative of orb. frequency, (Hz s“^) .. 5.97(5) x 10“^^ 

2nd derivative of orb. frequency, (Hz s“-) . 1.77(4) x 10“^^ 

3rd derivative of orb. frequency, (Hz s“^) . -1.06(2) x 10“^^ 

4th derivative of orb. frequency, (Hz s“^) . -5.6(2) x 10“^‘ 

5th derivative of orb. frequency, (Hz s'^) . 1.41(5) x 10“"^^ 

6th derivative of orb. frequency, (Hz s“^) . 1.15(6) x 10“^^ 

Epoch of ascending node, Tasc (MJD) . 55791.9182085(3) 

Projected semimajor axis, X (It-s) . 0.611656(4) 

Derivative of proj. semimajor axis, |jc| (It-s s“^) < 6 x 10“^“^ 

1st Laplace-Lagrange parameter, ei . 1(1) x 10“^ 

2nd Laplace-Lagrange parameter, €2 . —21(1) x 10“^ 


Derived Parameters 

RMS timing residual (fis) . 6.7 

Spin period, P (ms). 2.8842267415473(2) 

1st spin period derivative, P (s s*^). 1.4102(6) x 10“^® 

2nd spin period derivative, |P| (s“^) . <2x10“^* 

Orbital period, Porb (d). 0.1930984018(3) 

1 st derivative of orbital period, Porb (ss-‘) -1.66(1) X 

Eccentricity, e . 2.1(1) X 10^“* 


Gamma-ray Spectral Parameters 

Photon index, r . 1.2 it 0.3 

Cutoff energy, Ec (GeV) . 3.7 it 1.4 

Photon flux^ P ’100 (10“* photons cm“^ s**). 2.0it 0.2 

Energy flux^ Gioo (10^** erg cm“^ s“*). 2.8it0.1 


Note. — Numbers in parentheses are statistical Icr uncertainties in the 
last digits. The JPL DE405 solar system ephemeris has been used and times 
refer to TDB. 

“ Measured over the energy range from 100 MeV to 100 GeV. 


refined the pulse profile template. In the subsequent MC run, 
we then also let the sky location vary and determined the best- 
ht values for u, from which a further refined pulse prohle tem¬ 
plate was derived for use in the next MC run. Subsequently, 
we further extended the timing solution to also include addi¬ 
tional effects, such as proper motion and orbital eccentricity. 
We also tested for a variation in the projected semimajor axis, 
X, and a second spin-frequency derivative, /, but in both cases 
no signihcant measurement was made and so we can only give 
upper limits. 

Table presents the hnal timing solution, obtained after 
several iterations of pulsar-model and pulse-prohle-template 
rehnements. For the best estimate and uncertainty of each 
parameter in the timing model we report the mean value and 
standard deviation of the one-dimensional marginalized pos¬ 


terior distribution, as determined by the MultiNest algo¬ 
rithm. This timing model includes the polynomial expansion 
of the oscillating orbital frequency as given in Equation Q 
truncated at the sixth order (K = 6), which gave a maximum 
log£ value of about 580. 

We also repeated the timing procedure for different mod¬ 
els with polynomial expansions truncated at orders K = 5 and 
K =1. To determine which of the three models better fits 
the dat a, we then invo ked the Bayesian Information Criterion 
(BIC; |Schwarz|1978| l as a comparative tool. The BIC value 
is based on the number of free parameters to be estimated, 
the number of data points and the maximized value of log£. 
The smaller BIC value is preferred, since the BIC becomes 
increasingly large when either the data is poorly ht or when 
the number of free parameters increases and the data is over- 
ht. Thus, the BIC is minimized for the simplest model that 
sufficiently hts the data. In our case, the K = 1 model gave 
about the same log£ value as for K = 6, but given the extra 
parameter, the simpler K = 6 model is preferred according to 
the BIC. On the other hand, the K = 5 model has one param¬ 
eter less compared to K = 6, but gave a signihcantly lower 
log£ value, so that the BIC again preferred the K = 6 model. 
Therefore, we considered the model with K = 6 as shown in 
Table [T] as the optimal representation of the data. However, 
we also caution that outside the time span of this data, the de¬ 
rived timing solution, specihcally the polynomial expansion 
for the orbital frequency, likely has little predictive power. 


2.3. Timing Results 


In Figure [Tl the inferred parameters of the hnal timing so¬ 
lution of Tame [T] were used to generate the integrated pulse 
prohle (panel B), the diagram of rotational phase versus time 
(panel C), and the residuals over time (panel D). The recti- 
linearity of the phase tracks over time and the smallness of 
the residuals clearly certify the phase-coherence of the timing 
solution, accounting for the correct number of pulsar turns 
during the observational time span. Merely to illustrate the 
quality of the timing solution, in Figure (panel D) we also 
show the phase residuals that have an rms about 1.5% in turn, 
which translates into an rms timing accuracy of about 7 ^s. 

While the sky position of the hnal timing solution is 
compatible with the optical counterpart location, we were 
also able to measure signihcant values for the proper mo¬ 
tion of the system relative to the solar system barycen- 
ter {p,a,gLs, given in Table [T]). This amounts to a total 
transverse proper motion of p, = (/i^cos^(5-l-/r|)^/^ « (31 ± 
10) mas yr“'. Combined with the radio-DM distance es¬ 
timate of d = 0.45 kpc, we derive a transverse velocity of 
Vt = dgL, K, 30kms“*. This transverse movement contributes 


to the observe d spin parameter s P and P due to a changing 


Doppler shift (Shklovskii 1970 1 via; P = P\ni+Pv, /(dc). The 


latter term, representing the Shklovskii effect, is about 22% 
of the observed P, so that the intrinsic spin-period deriva¬ 
tive is estimated as Pim « 1.1 x 10“^° s s“^ Hence, we de¬ 
rive a Shklovskii-corrected pulsar spin-down luminosity of 
E = 47r^/Pint/F^ ~ 1-8 X 10^"^erg s“', where /= 10"^^gcm^ is 
an assumed hducial neutron star moment of inertia. 

The timing solution also provides constraints on the com¬ 
panion mass Me through the pulsar mass function by combin- 








































Gamma-ray Timing of PSR J2339-0533 


5 


ing the measurements of x and Porb as 




Ml sin^r 

= (Me+Mp)2=^ 

= (6.58942 ± 0.00012) X 


10 -^ Mo, 


(3) 


where G denotes the gravitational constant, Mp labels the pul¬ 
sar mass, and l is the inclination angle. When further com¬ 
bined with the optical observations of the companion mass 
estimates for bo th co mponents are possible, as will be eluci¬ 
dated in Section im 

We also measured an orbital eccentricity of 
e = 2.1(1) X 10“^. Given the short 4.6 hr orbital period 
of PSR J2339-0533, this value of residual eccentricity is 
higher than pre dicted by the convective-fluctuation theory of 
|Phinney| ( |1992| l. This may indicate additional evidence for 
convect ive f luctuations in the companion, as we discuss in 
Section [T4| While these ar e perhaps more complicated than 
assumed by|Phinney|(|1992|l, they could have inhibited perfect 
circularization of the binary orbit of PSR J2339-0533. 

Most remarkably, the timing solution unveils the extreme 
orbital-period variations of PSR J2339-0533 over the six 
years of LAT data. The two most extreme values of 
the orbital-period derivative attained during this interval are 
Porb = -5.8 X 10“'° s s-' and Porb = 2.7 x 10“'° s s”'. Panel E 
in Figure displays the dramatic evolution of the fractional 
changes APorb /Porb based on the ephemeris of Table A 
significant cyclic modulation of the orbital period is revealed, 
where the observational interval appears to cover about one 
and a half cycles. The possible causes of the Porb -modulation 
are discussed in detail in Section |3]below. 


3. DISCUSSION 
3.1. Component Masses 

By combining the opti cal radial velocity meas urements of 
the irradiated companion ( |Romani & Shaw 201 l| l and the pul¬ 
sar ephemeris, further system parameters can be constrained, 
such as the mass of the pulsar Mp and of the companion M^- 
In forthcoming work that will exploit the extremely accurate 
pulsar mass function when modeling the optical lightcurve, 
more precise limits on the system parameters will be p ossi¬ 
ble. However, here w e employ the previous estimates by |Ro-| 
|mani & Shaw| ( [ 20 TT] ) in order to obtain hducial values that are 
sufficient for the purpose of the subsequent discussion and 
conclusion as will be shown below. 


From radial velocity measurements of the companion, Ro 
mani & Shaw ( |2011| l observed a semi-amplitude of /fobs « 
270kms“F While these measurements track the apparent 
center of light, the radial-velocity amplitude o f the compan¬ 
ion’s center of mass, K 2 , can be larger (e.g., |van Kerkwijk| 
et al. 11201 l|l. The 'W-co rrection" for this effect inferred by 
Romani & Shaw |2011 1 gave K 2 « 1.18/fobs ~ 320kms“*, 
which we adopt here. The amplitude of the optical mod¬ 
ulation mainly depen ds on the inclination angle, which |Ro-| 
mani & Shaw| ( 20lT| estimated as i « 60°. As shown in Fig¬ 
ured this gives rise to the following estimates; Mp = 1.48 Mq 
ancUMc = 0.32Mq, that we assume as fiducial values for 
the rest of this paper. The resulting mass ratio would be 
^ = Mp/Mc = 4.61. 

For these mass estimates, the orbital separation is 
a « 1. 71 Rq and the co mpanion Roche lobe radius accord¬ 
ing to |Eg^eton| (|1983|l is: Ri « O. 88 R 0 . The actual ra¬ 
dius of the companion star, R^, is difficult to obtain, but can 
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Fig. 2.— Mass-mass diagram for the PSR J2339-0533 binary system. The 
non-shaded region is excluded by the pulsar mass function from the gamma- 
ray timing solution since sint < 1 (lower exclusion region) and from the 
companion mass function based on the optical radial-velocity measurements 
IRomani & Shaw]|20111 since K 2 > Xobs (upper exclusion region). Over 
the shaded (permitted) region, the color co de s hows how the fraction of the 
companion’s luminosity AL/L of Equation required by the gravitational 
quadrupole coupling model varies with the assumed component masses. The 
black dot indicates the fiducial values for the c omponent masses assum ed 
here, based on the optical lightcurve modeling by|Romani & Sha^l2011). 



also be es timated from the optical light-curve modeling. In 
doing so, [Romani & Shaw| ( 201 l| l found that the compan¬ 
ion is close to Roche lobe filling, with a filling factor of 
Rc/Rl ~ 0.9. Further support for this picture is provided by 
the observed radio eclipses during a large fraction of the or¬ 
bit ( [Ray et ar]|2014| 2015| l. Thus, this implies a companion 
radius of Rc « 0.79R©, which we also adopt as fiducial value 
for the remainder of this discus sion. _ 

The optical spectroscopy by [Romani & Shawj (pOTT) also 
shows that the companion’s spectrum resembles that of a late- 
type star. However, the above estimates suggest that the stellar 
companion is likely less dense than a main-sequence star of 
the same mass, which has also been found fo r other redback 
systems ( [Archibald et al.|2013[[Li et al.|2014| l. 

3.2. Energetics and Irradiation 

From the measured FAT energy flux Gioo (given in Ta- 
ble[T]), we can obtain an estimate of the gamma-ray luminosity 
= Aird^Gioa « 6.7 x 10^^ erg s“', assuming no geometric 
beaming correction. This leads to a gamma-ray conversion 
efficiency of = L^/E » 4%, not atypical compared to the 
MSPs studied in [Abdo et al.[ ( [2013| l. 

We can al so estimate the effect of p ulsar irradiation on the 
companion. [Romani & Shaw[ ([201 1[) inferred a companion- 
backside ternperature of approximately 2800 K, which sug¬ 
gests an intrinsic luminosity of L « 1.3 X 10^^ergs ' 
(0.034L©) based on the above estimated radius. The irradia¬ 
tion efficiency is often used to describe the fraction of the 
pulsar spin-down energy absorbed and converted into optical 
emission. Assuming the contribution of an isotropic pulsar 
wind to the companion’s optical luminosity via re-radiation 
is Lin- = 7yirrFR//(4a^) « 9.7 X lO^^pin erg s“*, we find an ef¬ 
ficiency of Pin = L/Lin « 14%. This is within the range of 
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typical ir radiation efficiencies as found for the systems inves¬ 
tigated in Breton et al.|(2013]l. 


3.3. Possible Causes for the Orbital Period Variations 
Generally, an observed rate of change in the binary orbital 
period, Pp rb, may result from various physical effects a s dis¬ 
cussed in ( |Doroshenko et al.|2001HLazaridis et al.|201 1[ ). The 
possible causes, most relevant regarding the time scale of our 
observations, can be summarized as 

(4) 

including energy loss due to gravitational-wave emission 
changes due to Doppler shifts (P°b)’ m^ss loss from 
the system (P^b)’ gravitational quadmpole moment cou¬ 
pling (^oa)- , 

The term P^* in Equation (Wl is the contribution due to 
gravitational-w ave emission. P^ot the case of circular orbits 
this is given by (|Peters|1964|l, 


pGW _ 
^oi-b - 


1927r 


27r 

Porb 


GM, 

r-i 


5/3 


q 

(^-b 1)1/3 


(5) 


where c denotes the speed of light. The resulting value for 
P™ ~ -1.3 X 10“'^ s s“i is about three orders of magnitude 
smaller than the measured values of Porb, and hence unlikely 
to be the primary cause. 

The term P^^ in Equation ([4 represents the combined effect 
of proper motion (|Shklovskn |T970| and Doppler shifts due 
to a changing distance to the binary system, e.g., from the 
Galactic acceleration or due to a massive third body; 

pD _ pShk I pGal . pace 

^orb “ “r 4 „rK 4 „i-K 


'^orb 


orb 


orb 


(6) 

Using the measured proper-motion values given in Table 
and the radio-based DM distance of d = 0.45 kpc, one obtains 


pShk _ 
^orb - 


(pi cotf 5+ til)d 


Port 1.8 X 10 “'‘*ss“i. 


(7) 


Contributions from the Gala ctic acceleration, P2-h ^ are typ¬ 
ically of similar magnitude (Lazaridis et al. 2009 1 . An ac¬ 
celeration of the binary system caused by a massive third 
body is also highly unlikely. If that was the case, the 
pulsar spin period and period derivative would be affected 
in the same way. Erom the measured values in Table 
we can estimate the maximum contribution of this effect 
by assuming the apparent spin-down is entirely due to ac¬ 


celeration; (P^b/Porb) = (P/P) — 5 X 10 ^®s \ implying 

^-14 c c -1 


Pgjb — 10 s s *. Thus, overall the term P^ 


is found to be 

more than four orders of magnitude smaller than the observed 
value for PSR J2339-0533. 

The term P^^ in Equation (|4 1 results from mass loss of the 
binary system. The spin-downTuminosity of the pulsar irradi¬ 
ating the companion can drive mass loss at a rate Me from the 
companion through an evaporative wind. Taking the orbit to 
be circu lar and assum ing no mass is lost from the pulsar, one 
obtains Peans|1924l l, 


^orb ' 


= -2—Porb, 


(8) 


where M = M^+ Mp denotes the total mass of the system. 
Assuming that all of the radiative energy from the pul¬ 
sar intercepted by the companion (for a presumed isotropic 


pulsar wind) is converted into mass loss from the system, 
this would imply a rate Me — 2.0 x lO“*M 0 yr“' ( Stevens 
et al. 1992| l. Inserting this value into Equation ([^ gives 
P^i, ~ -1.2 X 10“" s s“', which is more than an order of mag¬ 
nitude smaller in number than the measured Porb- More real¬ 
istically, adopting the 14% efficiency for this process as es- 
one gets Me — 2.8 x 10 “®MQyr“' and 
*. This therefore appears unlikely to 


3.2 


timated in Section 
P"b--l-7x lO-'^ss 
be the primary cause for the orbital variation observed. 

The possible contributions to Porb considered so far seem 
unlikely to be able to account for the large Porb-variation 
measured. In addition, these processes generate monotonic 
changes, which cannot explain the alternating increase and 
decrease of Porb observed. Hence this must be caused by the 
last term, P^j^, resulting from cyclic changes of the compan¬ 
ion’s gravitational quadrupole moment. In the following, we 
scrutinize the plausibility of this explanation, confronting a 
specific theoretical model with the observational data. 

3.4. Gravitational Quadrupole Coupling 

To explain the observed orbital-period modulation, the only 
plausible cause of those discussed above that remains is a 
changing gravitational quadrupole mo ment of the compan¬ 
ion star ([Matese 


& Whitmirel |1983|l. Such gravitational 
quadrupole coupling (GQC) li kely results from magnetic ac¬ 
tivity in the stell ar companion ([Applegate & Patterson|1987[ 

I Applegat^l 1 992| l . [Applegate & Shaham| ( |1994| l successfully 
applied the model by [Applegatej ( |1992l l to the black widow 
pulsar binary PSR B19574-20, gr avitationally coupling the 
orbital-period changes reported in ( [Arzoumanian et al.|1994 1 
to changes in the quadrupole moment of the companion star. 
Similar Porb variations on a time scale comparable to that 
found for PSR J2339-0533 were seen in another black widow 
system, PSR J2051-0827, and were also attribu ted to GQC 
( jPoroshenko et aL]|2001[ Lazaridis et al.j^llj ). The tran¬ 
sitional redback system PSR 110234-0038, which recently 
changed its state from MSP back to LMXB (jStappers et al.| 
2014|l, also displayed orbital-period changes at a comparable 


level ( [Archibald et al.|2013} . Eurther fitting into this picture, 
notably analog orbital-period changes were found in long 


term monitoring of LMXB systems (Wolff et al. 2009 Pa- 
truno et al.[|2012|l, likely also caused by GQC arising from 
cyclic magnetic activity. 

In Applegate’s GQC model the companion star is mag¬ 
netically active with a convective envelope and the pulsar is 
treated as a point mass moving in the gravitational field of 
the companion in a circular orbit. When the companion’s 
gravitational quadrupole moment (due to tidal and centrifugal 
deformations) varies with time, this can cause a variable or¬ 
bital motion at constant orbital angular momentum. When the 
companion becomes more oblate, its quadrupole moment in¬ 
creases, to balance gravitation the centripetal acceleration on 
the pulsar must increase, so that Porb must decrease, and vice 
versa. Therefore, in principle any mechanism that can modu¬ 
late the quadmpole moment of the companion also modulates 
the binary orbital period. _ 

The specific mechanism suggested by [Applegat^ \\992\ 
considers cyclic, solar-like magnetic activity m the convec¬ 
tive zone of the non-degenerate companion star. The mag¬ 
netic field within the companion is assumed to generate a 
torque, which cyclically exchanges angular momentum be¬ 
tween its outer convective layers and the inner part. Causing 
cyclic spin-up and spin-down of the companion’s outer layers. 
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this leads to an oscillation of the oblateness and the gravita¬ 
tional quadrupole moment of the companion. When the outer 
layers spin up, the star becomes more oblate, the quadrupole 
moment increases and Porb must decrease at fixed total orbital 
angular momentum. On the other hand, if the outer layers spin 
d own, the quadrupole moment decreases, and Porb increases. 
Applegate & Shaham ( 1994|l further hypothesized that tidal 


dissipation supplies the energy flows driving the convection 
in the rotating companion. The resulting magnetic activity 
and the mass loss due to the wind driven by pulsar irradiation 
contribute to a spin torque that holds the stellar companion 
slightly out of synchronous rotation giving rise to tidal dissi- 
p ation and thus heats the companion internally. 

Lanza et al. ( 1998[) proposed an extension of Applegate’s 


model. Instead of th e internal angular momentum redistri¬ 
bution considered by [Applegat^ \ \992\ , they argued that a 
variation of the azimuthal magrietic field can also produce 
quadtTipole moment changes. They invoked cyclic exchanges 
between kinetic and magnetic energy within the convective 
zone during the magnetic cycle to cause the modified distri¬ 
bution of angular momentum within the star. 

To better estimate the duration of the magnetic activity 
cycle coupling to the observed orbital-period modulation of 
PSR J2339-0533, we compared the timing solution to the 
following simple model. We described the change in orbital 
frequency A/orb(f) around its nominal value (according to 
Porb of Table [TJ by a term linearly varying with time t (mea¬ 
sured as difference from Pasc in Table [Til, and an oscillatory 
term with modulation period Pmod = l/^'mod and derivative 
Pmod — ^mod/^mod’ 


A/orb(f) = /mb.* f + A sin [l-nv^odt + 7ri>modf^ + ■/] , (9) 

where A and ip are the Poib-modulation amplitude and phase, 
respectively. Using chi-square minimization, the best-fit re¬ 
sults we found for this modulation model along with their for¬ 
mal Icr statistical uncertainties are given in Table [^ In partic¬ 
ular, we estimated a modulation period Pmod = AA yr and a 
modulation amplitude APc,rb/(27r) = APoib/Porb = 2.3 x 10“^. 
From the inferred value for /orb * we derived the residual 
orbital-period derivative Porb,* given in Table However, we 
caution that substantial systematic errors are expected, given 
that the reduced chi-square value of the fit was a few times 
greater than unity and secondary minima existed. Nonethe¬ 
less, as illustrated in panel E of Figure [T] this model follow¬ 
ing Equation (j^ allows one to qualitativ^y reproduce the ob¬ 
served orbital-period variation. 

These results are to be compared to the GQC model, where 
a change in orbital period by APorb is related to a change AQ 
of the companion’s quadrupole moment Q via (|Applegate &| 
|Patterson|1987| l, 


APo,b _ n ^6 

Porb Afc 


( 10 ) 


While Q likely has a complex dependency on how mass, 
magnetic fields and rotational angular velocity are distributed 
within the star, it is dominated by the mass distribution in the 
outer layers of the companion where the centrifugal acceler¬ 
ation is largest. To estimate the resulting change AQ from 
the transfer of ang ular momentum to the outer part of the star, 
Applegate ( 1992|l considered a thin shell of mass Ms and ra- 


dius Pc whose angular velocity U will change by Ail, giving 


TABLE 2 

Estimated Orbital-period Modulation Model Parameters 


Parameter 

Value 

Modulation amplitude, A (s“^). 

Modulation period, Pmod (yO. 

Modulation period derivative, Pmod (s s“^) .... 

Modulation phase, ijj (rad) . 

Residual orbital-period derivative, Porb.* (s 

8.48(2) X 10^'* 
4.15(1) 

-0.87(1) 

5.64(1) 

-1.57(2) X 10-“ 

Note. — Numbers in parentheses are only statistical, formal Icr uncer¬ 
tainties in the last digits. 

([Applegate & Shahamj 1994[|, 


M, Ail GM, APo,b 

(11) 

Me il 2Rl Rl 4^2 Roeb 

Assuming corotation of the tidally locked companion, we ob¬ 
tain for the PSR J2339-0533 system. 

Ail -i Me 

— =9.7x 10-’ 
il Ms 

(12) 


|Applegate| \\992\ found observational results to be typically 
ht well for a shell of mass Ms « O.lMc. Adopting this value 
here implies Ail/il « 10“^. In Applegate’s model, the mag¬ 
netic field generation and angular velocity variation operates 
cyclically with period Pmod (the cycle of the orbital-period 
modulation). Whi le no further deta i ls of this activity cy¬ 
cle we re specified, [Applegat^ ( |1992| l; [Applegate & Shahamj 
(|1994|) assumed that the energy needed for the transfer of an¬ 
gularmomentum would require an associated change in the 
star’s luminosity of 


AL = 


TT G Mp Ail AP 


orb 


3 Pmod Pc ^ Porb 


(13) 


Using the above estimated value for Pmod, the associated lumi¬ 
nosity change is AL « 4.1 x 10^^ ergs“*. This is only about 
0.03% of the intrinsi c lum inosity of the companion, which we 
estimated in Section [T2] Therefore, the companion should be 
easily capable of providing the required energy to power the 
orbital-period variations. Even if the true component masses 
were slightly different from those assumed, this conclusion 
remains robust. This is illustrated in Figure S where the 
required fractional luminosity is shown to varyby less than 
an order of magnitude over the permitted region in the mass- 
mass diagram. 

For other classes o f close binarie s, suc h as the RS C anum 
Venaticorum systems, |Lanza|(|2005|l (and |Lanza|2006|l found 
a discrepancy of Applegate’s GCQ model: The required sur¬ 
face angular velocity variations led to an energy dissipation 
rate in the turbulent convection zone of the secondary star ex¬ 
ceeding the stellar luminosity. For PSR J2339-0533, with 
lower APorb/Porb and Ail/il, this is not the case. Applying 
Equation (26) of|Lanza|p006|l, we infer a ratio between the 
dissipated power and the stellar luminosity of only a few per¬ 
cent, further supporting Applegate’s GCQ mechanism. How¬ 
ever, if the true companion radius was much smaller than the 
one we assumed above, this energy balance could eventually 
become questionable. But this would also imply that the com¬ 
panion was actually not near filling its R oche lobe, which 
is currently not supported observationally ([Romani & Shaw 

[20TT] |. 
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Moreover, on the timescale of the observational data, the 
component masses and the orbital angular momentum are 
approximately constant as presumed by the GQC model. 
This implies that the observed orbital-period variation must 
also lead to a change in the orbital separation: dja = 


2f’oib/^orb- Taking the largest Porb we observed, this would 
give i ~ 2 X 10“'^ It-s s“'. This is still smaller than the mea¬ 
sured upper limit on x provided in Tablej^and thus this aspect 
of the model is also in line with our timing measurements. 

The above analysis of the orbital-period evolution (Table|^ 
indicates that aside from the modulation, a long-term change 
of forb,* ~ 10“" s s“' might still remain. If this residual 
orbital-period derivative was mostly due to mass loss, the next 
largest contribution as estimated in Section 3.3 then the rate 
would be Me ~ 2 X 10~‘^MQyr~ ‘. To match this mass-loss rate 
the estimates made in Section 3.3 would require only little 
modification of the assumed irradiation efficiency or a slight 
deviation from the pulsar’s presumed isotropic emission. 

Overall, we conclude that the GQC theory offers a compat¬ 
ible explanation of the orbital-period variation observed for 


PSR J2339-0533. 


4. CONCLUSIONS 

Using the available Fermi-LAI data spanning more than six 
years, we carried out a precision gamma-ray timing analy¬ 
sis of the redback-type pulsar binary PSR J2339-0533. Most 
notably, the results revealed a long-term modulation of the 
4.6 hr binary orbital period. We found that this observed phe¬ 
nomenon can be explained by variations of the gravitational 
quadrupole moment of the co mpanion, through a mechnism 
proposed by |Applegate| ( |199^ . 

Since Applegate’s model offers a compatible explanation, it 
implicates that the companion star must have a sizeable outer 
convective zone, where cyclic magnetic activity causes the 
quadrupole moment changes leading to the observed orbital- 
period modulation. Additional evidence for such convective 
fluctuations is provided by our measured residual orbital ec¬ 
centricity, found to b e higher than th eoretically predicted for 
such a tight binary (|Phinney||1992|l. The strong irradiation 
and possibly also the tidal interaction with the close-by pulsar 
may drastically affect the internal structure of the companion 
making it similar to that of a main-sequence low-mass star 
with a convective envelope. While i ts optical spectrum is in 
fact consistent with a late-type star ( |Romani & Shaw|2011| l, 
the companion however appears to be less dense than a maiii- 
sequence star of the same mass. 

Combined with the pulsar ephemeris, future optical obser¬ 
vations of the companion can help to more tightly constrain 
the system parameters, e.g., their component masses. As an¬ 
other interesting prospect, such measurements might also be 
able to provide evidence for a variation in accordance with 


stellar activity. This might also allow further tests of the GQC 
theory, given the high timing precision with which we can 
measure the orbital-period variation. Ultimately, one might 
hope to identify the type of stellar magnetic dynamo in ac¬ 
tion or the physical origin of the differential rotation within 
the companion’s convective layers, currently observationally 
difficult to access otherwise. 

Surprisingly many black widow and redback pulsar binaries 
have been discovered in targeted radio searches of uniden¬ 
tified LAT sources. However, the fact that many of those 
systems also show significant radio eclipses and variations in 
their orbital parameters often makes it still extremely difficult 
to also detect the gamma-ray pulsations. That is because the 
short-term radio timing solutions cannot be extended imme¬ 
diately backwards to cover the full LAT data time span since 
FenriVs launch. As demonstrated for PSR J2339-0533, the 
method presented here can help address this problem. Besides 
revealing the gamma-ray pulsations, by making full use of the 
great potential of the exquisite LAT data, this also opens the 
unique possibility for long-term monitoring and insights into 
the complex interplay of effects influencing the orbital evolu¬ 
tion of such irradiating pulsar binary systems. 
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